#ifndef CCATACLISM_ACRETION_H
#define CCATACLISM_ACRETION_H

#include <list>
#include <algorithm>
#include <limits>
using namespace std;

//#define GRAVITY_CONSTANT 0.0000000000667428
#define GRAVITY_CONSTANT 6.67428e-11

template <typename _GAO> class CEmbryo {
public:
	typedef typename _GAO::gravity_t			gravity_t;
	typedef typename _GAO::velocity_t			velocity_t;
	typedef typename _GAO::container_t			container_t;
	typedef typename _GAO::pattern_t			pattern_t;
	typedef typename _GAO::label_t				label_t;
	typedef typename _GAO::planetesimal_t		planetesimal_t;
	typedef typename _GAO::particle_t			particle_t;

public:
	velocity_t		_velocity;
	container_t		_container;
	pattern_t		_representative;
	long double		_distance;
	long double		_covariance;

	CEmbryo (particle_t*& _p, pattern_t& _q) {
		_velocity = _p->_velocity;
		_representative = _p->_mean;
		_container = _p->_container;
		_distance = (_representative == _q) ? 0 : _representative - _q;
		_covariance = 0.0;
		container_t::iterator it = _container.begin();
		for (; it != _container.end(); it++)
			_covariance += _representative - *it;
		_covariance /= _container.size();
	}

	void update_values (const pattern_t& _query){
		pattern_t _rep(_representative.size());
		container_t::iterator it = _container.begin();
		for (; it != _container.end(); it++){
			_rep += *it;
		}
		_rep /= (double)_container.size();		
		_distance = (_rep == _query) ? 0 : _rep - _query;
		_representative = _rep;

		_covariance = 0.0;
		it = _container.begin();
		for (; it != _container.end(); it++)
			_covariance += _representative - *it;
		_covariance /= (double)_container.size();

	}

	bool operator < (const CEmbryo& _e){
		return _covariance < _e._covariance;
	}


	long double operator () (const CEmbryo& _e, double _g){
		long double _dividend = pow(_representative - _e._representative,2);
		if (_dividend == 0) return 0;
		return  _g * (_container.size() * _e._container.size()) 
				/ _dividend;
	}

	~CEmbryo (void){
	}

	template <typename __GAO> 
	friend ostream& operator<< (ostream& out, const CEmbryo<__GAO>& e){
		out << "Representativo: " << e._representative << endl;
		vector<string> _labels;
		vector<int> _counters;
		container_t::iterator it = e._container.begin();
		for (; it != e._container.end(); it++){
			bool exist = false;
			int k = 0;
			for (; k < _labels.size(); k++)
				if (_labels[k] == it->_label) {exist =  true; break;}

			if (exist) _counters[k] ++; 
			else {
				_labels.push_back(it->_label);
				_counters.push_back(1);
			}
		}
		for (int k = 0; k < _labels.size(); k++)
			cout << _labels[k] << " : " << _counters[k] << endl;
	}
};


template <typename _GAO> class  CCataclismAcretion {
public:	
	typedef typename _GAO::value_pattern_t		value_pattern_t; 
	typedef typename _GAO::gravity_t			gravity_t;
	typedef typename _GAO::velocity_t			velocity_t;
	typedef typename _GAO::container_t			container_t;
	typedef typename _GAO::pattern_t			pattern_t;
	typedef typename _GAO::label_t				label_t;
	typedef typename _GAO::planetesimal_t		planetesimal_t;
	typedef typename _GAO::particle_t			particle_t;
	typedef CEmbryo<_GAO>						embryo_t;
	typedef list<embryo_t>						protoplanet_t;
	typedef list<container_t>					feed_zone_t;

public:
	protoplanet_t	_protoplanets;
	protoplanet_t	_nprotoplanets;
	feed_zone_t 	_q_zones;
	pattern_t		_query;


	void __export (){
		feed_zone_t::iterator fid = _q_zones.begin();
		
		for (; fid != _q_zones.end(); fid++){
			ofstream out ("plot.txt");
			ofstream lbls ("labels.txt");
			double s_angle = 360.00 / (double) fid->size();
			container_t::iterator pit = fid->begin();
			pattern_t _mean = pattern_t(pit->size());
			for (; pit != fid->end(); pit++){
				_mean += *pit;
			}
			_mean /= fid->size();

/*			
			pit = fid->begin();
			double _angle = 0.0;
			for (; pit != fid->end(); pit++){
				float x = (_mean - *pit) * cos(_angle);
				float y = (_mean - *pit) * sin(_angle);
				out << x * 10 << "\t" << y * 10 << " " << pit->_label << endl;
				_angle += s_angle;
			}
			*/

			pit = fid->begin();
			for (; pit != fid->end(); pit++){
				out << *pit - _mean << endl;
				lbls << pit->_label << endl;

			}
			out.close();
			lbls.close();
			system("@start plot.txt");
			system("pause");
		}
	}


	CCataclismAcretion (planetesimal_t& _planetesimals, const pattern_t& _q){
		_query = _q;
		planetesimal_t::iterator it = _planetesimals.begin(); 
		for (; it != _planetesimals.end(); it++){
			_protoplanets.push_back(embryo_t(*it,_query));
		}
		_protoplanets.sort();
	}

	CCataclismAcretion (planetesimal_t& _planetesimals, vector<pattern_t>& _queries){
		vector<pattern_t>::iterator qi = _queries.begin();
		for (; qi != _queries.end(); i++)
			_query += *qi;
		_query /= _queries.size(); 

		planetesimal_t::iterator it = _planetesimals.begin(); 
		for (; it != _planetesimals.end(); it++){
			_protoplanets.push_back(embryo_t(*it,_query));
		}
		_protoplanets.sort();
	}

	void operator () (int _t, int _z, double _b) {
		ofstream map ("map.txt");
		ofstream evol("evol.txt");
		int _time_period = 1;
		while (_time_period < _t + 1){
			double _gc =  GRAVITY_CONSTANT + pow((1 / _time_period),_b);

			cout << "Iteracion " << _time_period  << endl;
			double prom = 0.0;
			protoplanet_t::iterator pit = _protoplanets.begin(); 
			for (; pit != _protoplanets.end(); pit++){
				prom += pit->_covariance;
			}
			prom /= (double)_protoplanets.size();
			cout << "Promedio de evolucion : " << prom << endl;
			evol << prom << endl;
			_protoplanets.sort();

			pit = _protoplanets.begin(); 
			while (pit != _protoplanets.end()){
				__cross_over(pit,_gc);
				if (!pit->_container.empty())
					_nprotoplanets.push_back(*pit);
				protoplanet_t::iterator tpit = pit; pit++;
				_protoplanets.erase(tpit); 
			}
			_protoplanets = _nprotoplanets;
			_nprotoplanets.clear();
			_time_period ++;
		}
		evol.close();
		protoplanet_t::iterator  pit = _protoplanets.begin(); 
		for (; pit != _protoplanets.end(); pit++)
			swap (pit->_distance, pit->_covariance);
		_protoplanets.sort();

		__make_feed_zones(_z);
		// Imprimimos 
		feed_zone_t::iterator	fid = _q_zones.begin();
		for (; fid != _q_zones.end(); fid++){
			cout << "\nZone: " << endl;
			vector<string> _labels;
			vector<int> _counters;
			container_t::iterator it = fid->begin();
			for (; it != fid->end(); it++){
				bool exist = false;
				int k = 0;
				for (; k < _labels.size(); k++)
					if (_labels[k] == it->_label) {exist =  true; break;}

				if (exist) _counters[k] ++; 
				else {
					_labels.push_back(it->_label);
					_counters.push_back(1);
				}
			}
			for (int k = 0; k < _labels.size(); k++)
				cout << _labels[k] << " : " << _counters[k] << endl;
		}
	}
	
	void __cross_over (typename protoplanet_t::iterator& it, double _g){
		if(_protoplanets.size() == 1) return;
		protoplanet_t::iterator nit = _protoplanets.begin();
		protoplanet_t::iterator sit = _protoplanets.begin();
		long double _gravity = (*it)(*nit,_g); nit++;
		for (; nit != _protoplanets.end(); nit++){
			long double _ngravity = (*it)(*nit,_g);
			if (_ngravity > _gravity){
				sit = nit;
				_gravity = _ngravity;
			}
		}

		typename embryo_t::container_t min_container, max_container;
		typename embryo_t::container_t::iterator eit = it->_container.begin();

		for (; eit != it->_container.end(); eit++)
			sit->_container.push_back(*eit);
		sit->update_values(_query);

		eit = sit->_container.begin();
		for (; eit != sit->_container.end(); eit++)
			if (sit->_representative < *eit) max_container.push_back(*eit);
			else min_container.push_back(*eit);

		it->_container.clear(); 
		it->_container = min_container;
		it->update_values(_query);
//		_protoplanets.erase(sit);
/*
		// Datos exparcidos en el espacio tiempo 
		eit = max_container.begin();
		for (; eit != max_container.end(); eit++){
			nit = _protoplanets.begin();
			sit = _protoplanets.begin();
			double min_dist = *eit - nit->_representative;	nit++;
			for(; nit != _protoplanets.end(); nit++){
				double new_min = *eit - nit->_representative;
				if (new_min < min_dist) {
					sit = nit; min_dist = new_min;
				}
			}
			sit->_container.push_back(*eit);
			sit->update_values(_query);
		}
*/
		sit->_container.clear(); 
		sit->_container = max_container;
		sit->update_values(_query);
	}

	void __make_feed_zones (int _z){
		double lambda = (double)_protoplanets.size() / _z;
		protoplanet_t::iterator it = _protoplanets.begin();

		for (int i = 0; i < _z; i++){
			container_t _res_container;
			for (int k = 0; k < lambda; k++) {
				if (it == _protoplanets.end()) break;
				container_t::iterator pit = it->_container.begin();
				for (; pit != it->_container.end(); pit++)
					_res_container.push_back(*pit);
				it++;
			}
			_q_zones.push_back(_res_container);
		}
	}
};

#endif 
